Quantitative Big Imaging

Kevin Mader
20 March 2014

ETHZ: 227-0966-00L

Analysis of Complex Objects

Course Outline

  • 20th February - Introductory Lecture
  • 27th February - Filtering and Image Enhancement (A. Kaestner)
  • 6th March - Basic Segmentation, Discrete Binary Structures
  • 13th March - Advanced Segmentation
  • 20th March - Analyzing Single Objects
  • 27th March - Analyzing Complex Objects
  • 3rd April - Spatial Distribution
  • 10th April - Statistics and Reproducibility
  • 17th April - Dynamic Experiments
  • 8th May - Big Data
  • 15th May - Guest Lecture - Applications in Material Science
  • 22th May - Project Presentations

Literature / Useful References

  • Jean Claude, Morphometry with R
  • Online through ETHZ
  • Buy it
  • John C. Russ, “The Image Processing Handbook”,(Boca Raton, CRC Press)
  • Available online within domain ethz.ch (or proxy.ethz.ch / public VPN)

Previously on QBI ...

  • Image Enhancment
    • Highlighting the contrast of interest in images
    • Minimizing Noise
  • Segementation
    • Understanding value histograms
    • Dealing with multi-valued data
  • Automatic Methods
    • Hysteresis Method, K-Means Analysis
  • Regions of Interest
    • Contouring
  • Component Labeling
  • Single Shape Analysis

Outline

  • Motivation (Why and How?)
  • What are Distance Maps?
  • Skeletons
    • Tortuosity
  • What are thickness maps?
    • Thickness with Skeletons
  • Watershed Segmentation
    • Connected Objects
  • Curvature
    • Characteristic Shapes

Learning Objectives

Motivation (Why and How?)

  • How do we measure distances between many objects?
  • How can we extract topology of a structure?

  • How can we measure sizes in complicated objects?

  • How do we measure sizes relavant for diffusion or other local processes?

  • How do we identify seperate objects when they are connected?

  • How do we investigate surfaces in more detail and their shape?

  • How can we compare shape of complex objects when they grow?

    • Are there characteristic shape metrics?

What did we want in the first place

To simplify our data, but an ellipse model is too simple for many shapes

So while bounding box and ellipse-based models are useful for many object and cells, they do a very poor job with the sample below.

Single Cell

Why

  • We assume an entity consists of connected pixels (wrong)
  • We assume the objects are well modeled by an ellipse (also wrong)

What to do?

  • Is it 3 connected objects which should all be analzed seperately?
  • If we could divide it, we could then analyze each spart as an ellipse
  • Is it one network of objects and we want to know about the constrictions?
  • Is it a cell or organelle with docking sites for cell?
  • Neither extents nor anisotropy are very meaningful, we need a more specific metric which can characterize

Distance Maps: What are they

A map (or image) of distances. Each point in the map is the distance that point is from a given feature of interest (surface of an object, ROI, center of object, etc)

Simple Circles

Definition

If we start with an image as a collection of points divided into two categories

  • \( Im(x,y)= \) {Foreground, Background}
  • We can define a distance map operator (\( dist \)) that transforms the image into a distance map

\[ dist(\vec{x}) = \textrm{min}(||\vec{x}-\vec{y}|| \forall \vec{y} \in \textrm{Background}) \]

We will use Euclidean distance \( ||\vec{x}-\vec{y}|| \) for this class but there are other metrics which make sense when dealing with other types of data like Manhattan/City-block or weighted metrics.

Distance Maps: Precaution

The distance map is one of the crictical points where the resolution of the imaging system is important.

  • We measure distances computationally in pixels or voxels
  • but for them to have a meaning physically they must be converted
    • Isotropic imaging (1 \( \mu \) m x 1\( \mu \) m x 1 \( \mu \) m) is fine

Anisotropic

Ideally

  • as part of filtering, resample and convert to an isotropic scale.

If that is not possible

  • custom distance map algorithms which use the side-lengths of the voxels to calculate distance rather than assuming 1x1x1

Distance Maps

Simple Circles

We can create 2 distance maps

  1. Foreground \( \rightarrow \) Background

    • Information about the objects size and interior
  2. Background \( \rightarrow \) Foreground

    • Information about the distance / space between objects

Simple Circles

Simple Circles

Distance Map

One of the most useful components of the distance map is that it is relatively insensitive to small changes in connectivity.

  • Component Labeling would find radically different results for these two images
    • One has 4 small circles
    • One has 1 big blob

Simple Circles

Closer Circles

Circles Histogram

Distance Map of Cell

Foreground

Cell Distance Map

Background

Cell Distance Map

Skeletonization / Networks

For some structures like cellular materials and trabecular bone, we want a more detailed analysis than just thickness. We want to know

  • which structures are connected
  • how they are connected
  • express the network in a simple manner
    • quantify tortuosity
    • branching

Mesh Mask

Skeletonization

The first step is to take the distance transform the structure \[ I_d(x,y) = \textrm{dist}(I(x,y)) \] We can see in this image there are already local maxima that form a sort of backbone which closely maps to what we are interested in.

Distance Map

By using the laplacian filter as an approximate for the derivative operator which finds the values which high local gradients.

\[ \nabla I_{d}(x,y) = (\frac{\delta^2}{\delta x^2}+\frac{\delta^2}{\delta y^2})I_d \approx \underbrace{\begin{bmatrix} -1 & -1 & -1 \\ -1 & 8 & -1 \\ -1 & -1 & -1 \end{bmatrix}}_{\textrm{Laplacian Kernel}} \otimes I_d(x,y) \]

Skeleton

Creating the skeleton

We can locate the local maxima of the structure by setting a minimum surface distance \[ I_d(x,y)>MIN-SLOPE \] and combining it with a minimum slope value \[ \nabla I_d(x,y) > MIN-DIST \]

Thresholds

Harking back to our first lectures, this can be seen as a 2D threshold of the entire dataset.

  • We first make the dataset into a tuple

\[ \textrm{cImg}(x,y) = \langle \underbrace{I_d(x,y)}_1, \underbrace{\nabla I_d(x,y)}_2 \rangle \]

\[ \textrm{skelImage}(x,y) = \] \[ \begin{cases} 1, & \textrm{cImg}_1(x,y)\geq MIN-DIST \\ \textbf{ AND} & \textrm{cImg}_2(x,y)\geq MIN-SLOPE \\ 0, & \textrm{otherwise} \end{cases} \]

Skeleton: Different Thresholds

Skeleton

Skeleton

Skeleton: Pruning

Skeleton

Pruning

  • The structure is a overgrown
  • Stricter 'thresholds'

Skeleton: Junctions

With the skeleton which is ideally one voxel thick, we can characterize the junctions in the system by looking at the neighborhood of each voxel.

Skeleton

As with nearly every operation, the neighborhood we define is important and we can see that if we use a box neighborhood vs a cross neighborhood (4 vs 8 adjacent) we count diagonal connections differently

Skeleton

Skeleton: Network Assessment

Once we have classified the skeleton into terminals, segments, and junctions we can analyze the various components and assess metrics like connectivity, branching, and many other.

The easist way is to further process the image by

  • component labeling the different segments and joints
  • determining the average length of each segment
  • number of incoming segments at a given joint

Most of the other metrics can be simply counted

Cross Neighborhood

label x
Junction 720
Segment 332
Terminal 558

Box Neighborhood

label x
Junction 1321
Segment 783
Terminal 379

Skeleton: Tortuosity

One of the more interesting ones in material science is called tortuosity and it is defined as the ratio between the arc-length of a segment and the distance between its starting and ending points. \[ \tau = \frac{L}{C} \]

A high degree of tortuosity indicates that the network is convoluted and is important when estimating or predicting flow rates. Specifically

  • in geology it is an indication that diffusion and fluid transport will occur more slowly
  • in analytical chemistry it is utilized to perform size exclusion chromatography
  • in vascular tissue it can be a sign of pathology.

Thickness Map: What is it?

Thickness Map

Cell Distance Map

Cell Distance Map

Thickness Map

Cell Distance Map

Cell Distance Map

Thickness Map From Skeleton

Calculating the distance map by drawing a sphere at each point is very time consuming (\( O(n^3) \)).

  • The skeleton (last section) is very closely related to the thickness.
  • We found the local maxima in the image using the Laplace
  • We can thus grow the Spheres from these points instead of all
  • Start by instead of thresholding transforming the image to the distance at each point

\[ \textrm{thSkelImage}(x,y) = \] \[ \begin{cases} \textrm{cImg}_1(x,y) , & \textrm{cImg}_1(x,y)\geq MIN-DIST \\ & \& \textrm{ cImg}_2(x,y)\geq MIN-SLOPE \\ 0, & \textrm{otherwise} \end{cases} \]

Skeleton

From Skeleton vs All Points

   user  system elapsed 
  5.494   0.114   5.759 

Thickness from skeleton

   user  system elapsed 
 16.941   0.466  17.812 

Thickness from all points

Statistically Does it Matter

Thickness from all points

It depends

  • Small structures are lost
  • They might not have been very important or noisy anyways
  • Higher values are very similar
id Full Map Skeleton Map
Min. 0.75 1.75
1st Qu. 2.47 2.47
Median 2.66 2.66
Mean 2.67 2.74
3rd Qu. 3.01 3.01
Max. 3.20 3.20

How much can we cut down

   user  system elapsed 
  3.100   0.066   3.605 

Skeleton

Thickness Distributions Compared

id Full Map Skeleton Map Tiny Skeleton Map
Min. 0.75 1.75 2.02
1st Qu. 2.47 2.47 2.50
Median 2.66 2.66 2.66
Mean 2.67 2.74 2.77
3rd Qu. 3.01 3.01 3.01
Max. 3.20 3.20 3.20

Thickness in 3D Images

While the examples and demonstrations so far have been shown in 2D, the same exact technique can be applied to 3D data as well. For example for this liquid foam structure

Liquid Foam Plateau Border

  • The thickness can be calculated of the background (air) voxels in the same manner.

  • With care, this can be used as a proxy for bubble size distribution in systems where all the bubbles are connected to it is difficult to identify single ones.

Liquid Foam Thickness

Watershed

Watershed is a method for segmenting objects without using component labeling.

  • It utilizes the shape of structures to find objects
  • Imagine how rain would fall on a topographical map

Closer Circles

Curvature

“Computation of Surface Curvature from Range Images Using Geometrically Intrinsic Weights”*, T. Kurita and P. Boulanger, 1992.

http://en.wikipedia.org/wiki/Mean_curvature

Characteristic Shape

Curvature: Take Home Message

It provides another metric for characterizing complex shapes

  • Particularly useful for examining interfaces
    • Folds, saddles, and many other types of points are not characterized well be ellipsoids or thickness maps

Other Techniques

There are hundreds of other techniques which can be applied to these complicated structures, but they go beyond the scope of this course. Many of them are model-based which means they work well but only for particular types of samples or images. Of the more general techniques several which are easily testable inside of FIJI are

  • Directional Analysis = Looking at the orientation of different components using Fourier analysis (Analyze \( \rightarrow \) Directionality)
  • Tubeness / Surfaceness (Plugins \( \rightarrow \) Analyze \( \rightarrow \)) characterize binary images and the shape at each point similar to curvature but with a different underlying model
  • Fractal Dimensionality = A metric for assessing the structure as you scale up and down by examining various spatial relationships
    • Ma, D., Stoica, A. D., & Wang, X.-L. (2009). Power-law scaling and fractal nature of medium-range order in metallic glasses. Nature Materials, 8(1), 30–4. doi:10.1038/nmat2340
  • Two (or more) point correlation functions = Used in theoretical material science and physics to describe random materials and can be used to characterize distances, orientations, and organization in complex samples
    • Jiao, Y., Stillinger, F., & Torquato, S. (2007). Modeling heterogeneous materials via two-point correlation functions: Basic principles. Physical Review E, 76(3). doi:10.1103/PhysRevE.76.031110
    • Andrey, P., Kiêu, K., Kress, C., Lehmann, G., Tirichine, L., Liu, Z., … Debey, P. (2010). Statistical analysis of 3D images detects regular spatial distributions of centromeres and chromocenters in animal and plant nuclei. PLoS Computational Biology, 6(7), e1000853. doi:10.1371/journal.pcbi.1000853
    • Haghpanahi, M., & Miramini, S. (2008). Extraction of morphological parameters of tissue engineering scaffolds using two-point correlation function, 463–466. Retrieved from http://portal.acm.org/citation.cfm?id=1713360.1713456